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•  INTRODUCTION 


Chemical  vapor  deposition  (CVD)  is  one  of  the  most  important  techniques  for  the 
processing  of  structural  materials.  It  is  extremely  versatile  in  that  a  wide  range  of 
materials  and  coatings  is  done  by  CVD,  perhaps  more  than  by  any  other  technique.  The 
process  is  flexible  enough  to  accommodate  a  variety  of  sizes  and  shapes  of  objects  to  be 
processed  and  allows  continuous  change  of  composition  of  the  material  or  coating  [1].  In 
this  report  we  discuss  the  application  of  CVD  processes  to  ‘structural  materials’.  The 
classification  ‘structural  materials’  here  is  based  on  the  end  use  being  in  aerospace, 
mechanical  and  related  applications  (as  opposed  to  electronics  materials)  and  hence  covers 
coatings  and  high  temperature  aerospace  materials  such  as  composite  structures,  ceramic- 
matrix  composites  (CMC),  functionally-gradient  materials  (FGM),  etc.  Representative 
examples  are  SiC  (fiber,  coating),  Si3N4,  BN,  other  carbides  and  nitrides,  composites  such 
as  SiC-SiC,  carbon-SiC,  carbon-carbon,  and  carbon-BN  [1]. 

There  are  several  variations  of  CVD  in  the  preparation  of  structural  materials. 
These  include: 

(1)  CVD  coating  of  complex  arbitrary  shapes  and  internal  cavities  such  as  occur  in 
the  following  examples.  Hard  CVD  coatings  are  used  in  rocket  nozzles  and  gas 
turbine  engines  for  erosion  and  wear  resistance.  T1B2  and  TiC  coatings  are  used 
for  coating  titanium  blades  and  Si3N4  coatings  on  carbon  nozzles.  AI2O3,  TiN,  TiC 
and  TiB2  coatings  are  used  on  tool  bits  to  increase  the  cutting  life  of  tools  through 
improvements  in  hardnesis  and  wear  resistance.  CVD  metal  films  are  used  as 
protective  coatings. 

(2)  Conventional  CVD  on  substrates;  This  offers  a  proven  method  for  producing 
dense  crystalline  materials  for  a  variety  of  applications;  monolithic  ceramic 
materials  such  as  SiC,  TiC,  boron;  ceramic  materials  as  diffusion  barriers. 

(3)  CVD  of  fibers.  Examples  include  SiC,  Si3N4,  boron,  and  oxide  fibers.  These 
monofilament  fibers  are  used  primarily  for  reinforcing  other  materials  in  producing 
desirable  composites. 

(4)  CVD  of  preforms.  Here  the  CVD  process  is  used  to  infiltrate  a  preform  made 
from  one  of  the  above  fibers  to  form  a  composite.  CVD  inside  preforms  is 
generally  called  chemical  vapor  infiltration  (CVI).  Composites  such  as  SiC-SiC, 
C-C,  C-SiC,  and  C-BN  have  been  produced  by  this  technique. 

(5)  Plasma-based  CVD.  This  variation  allows  lower  processing  temperatures  that 
may  be  critical  in  some  applications  to  obtain  the  desired  microstructures.  While  it 
is  popular  in  semiconductor  processing,  plasma-based  CVD  has  been  recently 
considered  in  processing  of  structural  materials. 


Materials  processed  by  CVD  have  found  a  large  number  of  applications  in 
aerospace,  automotive,  mechanical,  and  related  industries.  In  addition  to  the  applications 
mentioned  above  other  applications  include  SiC  produced  by  CVD.  SiC  produced  by 
CVD  alone  commands  a  range  of  applications.  For  example,  CVD-grown  SiC  is  used  for 
high  temperature  engineering  applications  such  as  reaction  tubes,  furnace  components  and 
liners,  heating  elements,  and  refractory  ware.  CVD  SiC  is  a  good  candidate  for  materials 
that  must  withstand  severe  environments.  This  includes  transmissive  domes  in  high-speed 
missiles,  reflective  and  transmissive  optical  components,  UV  telescopes  and  instruments 
for  astronomy  and  chemically  inert  industrial  seals.  SiC  is  also  attractive  for  storage  media 
heads  and  disks  for  high  packing  density,  and  finally,  in  semiconductor  device  applications. 

Carbon-carbon  composites  are  used  in  rocket  nozzles,  nose  cones  for  re-entry 
vehicles,  heat  shields,  and  aircraft  brakes  [2].  SiC-SiC  and  C-SiC  find  many  high 
temperature  applications.  Zr02  composite  reinforced  with  TiN  whiskers  is  a  good 
candidate  for  advanced  heat  engines  due  to  a  good  match  of  thermal  properties.  In 
producing  some  composites,  improvements  in  reliability  and  performance  can  be  expected 
when  such  materials  are  designed  with  a  continuous  compositional,  and  therefore, 
functional  gradient.  In  summary,  application  for  CVD-grown  materials  is  rather  extensive 
and  it  is  not  an  overstatement  to  say  that  the  related  commerce  including  end  uses  (both 
structural  and  electronics)  can  be  measured  in  billions  of  dollars. 

Despite  the  impressive  range  of  applications  outlined  above,  a  fundamental 
understanding  of  CVD  of  structural  materials  is  lacking.  A  recent  National  Research 
Council  report  [3]  on  high  performance  composites,  prepared  by  a  distinguished  panel  for 
the  National  Materials  Advisory  Board,  states;  “A  basic  science  approach  to  the  CVD..., 
providing  fundamental  insights  into  decomposition  kinetics,  fiber-formation  mechanisms, 
and  surface  reactions,  can  result  in  a  better  optimization  of  processes  and  products. 
Although  it  may  not  be  inordinately  difficult  to  form  a  new  CVD  fiber  of  a  given 
composition  initially,  achievement  of  consistency  in  properties  may  require  understanding 
and  carefully  controlling  the  parameters  of  the  deposition  process.” 

The  above  conclusion  is  even  more  valid  when  the  composition  is  varied  to 
produce  FGMs.  In  addition,  process  scale-up  from  laboratory  to  large  scale 
manufacturing  is  difficult  because  of  the  complexity  of  the  CVD  process  which  is 
characterized  by  the  simultaneous  heat,  momentum  and  mass  transfers  along  with  gas 
phase  (homogeneous)  and  surface  (heterogeneous)  reactions  and  complex  flow-induced 
orientation  effects.  In  this  regard,  mathematical  modeling  of  the  transport  phenomena 
associated  with  the  process,  and  computer  simulation  is  valuable  in  many  respects, 
including; 

(1)  providing  an  understanding  of  the  process  mechanisms; 

(2)  developing  a  cause  and  effect  relationships  between  the  process  variables  (for 
which  there  are  knobs  on  the  reactor  panel)  vs.  process  figures-of-merit,  such  as 
growth  rate,  uniformity,  surface  material  composition,  etc.  (note  that  these  figures- 
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of-merit,  in  turn,  affect  the  structural  characteristics,  including  mechanical  and 
thermal  properties); 

(3)  process  optimization  to  meet  production  goals; 

(4)  equipment  design  to  meet  given  processing  objective;  and 

(5)  process  control. 

Although  models  and  computer  code  can  serve  as  a  cost-effective  design  tool,  to 
date  there  has  been  very  little  modeling  work  for  CVD  of  composites.  This  lack  of 
application  results  from  a  variety  of  issues  regarding  physical  models  and  computational 
techniques  that  must  be  resolved  to  allow  practical  application  of  state-of-the-art 
computational  tools  to  the  problem  of  CVD  formed  composites.  These  include: 

•  models  for  multi-component  mass  transfer  with  chemical  reactions  in  composites 
processing, 

•  surface  reactions  and  related  surface  phenomena, 

•  sound  models  for  thermal  diffusion, 

•  preform  infiltration, 

•  model  to  represent  the  random  pore  network  of  the  preforms, 

•  development  of  suitable  boundary  conditions, 

•  transient  analysis  to  study  growth  of  FGMs 

•  computational  techniques  appropriate  for  the  above  physical  processes  both  in  terms 
of  simulation  speed  and  accuracy,  and 

•  models  to  correlate  the  transport  phenomena  effects  to  the  microstructure,  mechanical, 
and  thermal  properties. 

The  Small  Business  Technology  Transfer  (STTR),  that  is  reported  herein, 
addresses  these  issues  as  it  develops  suitable  transport  and  reaction  models  for  CVD  and 
CVI  of  composites  within  the  framework  of  a  computational  simulation  tool.  This  work 
has  significant  relevance  to  Air  Force  missions  since  high  temperature  composites  are 
extensively  used  in  military  applications.  Affordability  and  reliability  of  materials  are  two 
ultimate  factors  that  determine  the  success  of  modem  weapon  systems.  In  general,  the 
required  materials  are  very  specialized,  low  volume,  and  come  with  stringent 
specifications.  Improved  processes  through  a  better  understanding  provided  by  modeling 
will  impact  on  the  material  quality  and  cost  of  processing. 
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•ANALYSIS 


CVD  Processes 

Modeling  of  the  CVD  processes  require  solution  of  the  appropriate  governing 
equations,  specification  of  appropriate  boundary  conditions,  specification  of  mathematical 
models  of  relevant  physical  process  and  an  efficient  and  accurate  solution  procedure. 
Those  developed  and  used  in  the  present  effort  are  discussed  in  this  Section. 


Governing  Equations 

At  pressures  and  temperatures  normally  used  in  CVD  of  structural  materials,  the 
mean  fi-ee  path  of  the  gas  molecules  is  much  smaller  than  the  reactor  dimensions  and  as 
such,  a  continuum  analysis  is  appropriate  to  study  the  transport  and  kinetics  inside  the 
reactor.  The  following  are  the  global  mass,  momenta,  energy  and  individual  species  mass 


conservation  equations,  respectively. 

^  +  V-pU  =  0 
dt 

(1) 

^^^-i-V-pUU  =  -VP  +  V-x  +  pg 
dt 

(2) 

^^^+V-phU  =  -V-q+  — +  0  +  Q.,, 
at  Dt 

(3) 

^+Vpa>,U  =  -V.j,+SR,,  (i  =  l,N) 

(4) 

where  p  is  the  total  density,  U  is  the  mass-averaged  gas  velocity,  P  is  the  pressure,  x  is  the 
viscous  stress  tensor  and  g  is  the  gravitational  vector,  h  is  the  static  enthalpy,  q  is  the 
total  energy  flux  relative  to  mass  average  velocity  (i.e.,  diffiisional  energy  flux)  and  <I>  is 
the  mean  flow  dissipation  rate.  Qext  represents  external  heat  sources  including  radiation,  ji 
is  the  diffusive  mass  flux  of  the  i*  species  and  cO;  =Pi  /p  is  the  mass  fraction  of  the  i* 

J 

species.  The  last  term  in  Eq.  (4),  ’  represents  the  i*  species  mass  production  or 

j=i 

loss  through  the  J  homogeneous  chemical  reactions. 

The  above  transport  equations  need  to  be  augmented  with  an  equation  of  state, 
viz.,  the  perfect  gas  law. 


4 


M 


(5) 


where  T  is  the  gas  temperature  and  R  is  the  universal  gas  constant.  Here  the  mixture 
molecular  weight,  M,  is  given  by 


M  = 


(6) 


where  N  is  the  number  of  gaseous  species  and  M;  is  the  molecular  weight  of  the  i'*’  species. 
Further  details  on  some  of  the  key  terms  in  the  above  equations  are  given  below. 

The  static  enthalpy  per  unit  mass,  h,  is  the  mixture  enthalpy  of  N  species. 


h  =  2;a>,hi(T) 

i=l 


(7) 


where  hj(T)  is  the  enthalpy  per  unit  mass  of  the  i*  species  (a  function  of  temperature  only). 
The  temperature,  T,  is  related  to  enthalpy,  h,  via  the  relationship 


i'=2:«>,k  +  CcjT’)dr 

■  i=1  L  J 


(8) 


where  h^,  is  heat  (enthalpy)  of  formation  of  the  i*  species  at  temperature  T^,  and  Cp.  (T  ) 
is  the  specific  heat  at  constant  pressure  for  the  i***  species. 

The  total  energy  flux,  q,  is  the  sum  of  the  conductive  (Fourier)  and  inter- 
diffusional  contributions,  qc  and  qa,  respectively  where 


q,  =-kVT 

(9) 

N 

qd=S>'i(T)j, 

(10) 

i=l 


Here  k  is  thermal  conductivity. 

The  viscous  stress  tensor,  x,  is  related  to  the  rate  of  strain  tensor,  D,  and  the 
molecular  viscosity,  p,  by  the  relationship 


T  =  p 


2D-|(V-U)5 


(11) 
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where 


(12) 


D  =  1(VU  +  VU^) 

where  the  superscript  t  refers  to  the  transpose.  The  dissipation,  O,  is  given  by 

0  =  2pD:D-|p(V-Uf  (13) 

For  the  results  reported  in  this  study,  the  Reynolds  numbers  were  sufficiently  low  enough 
that  the  effects  of  turbulence  were  negligible.  Thus  the  Reynolds  stress  terms  were 
neglected. 

Molecular  Properties 

For  the  cases  of  interest  in  this  study  the  gases  consist  of  mixtures  of  individual 
gaseous  species.  Therefore  it  was  necessary  to  calculate  multiple  species  values  of 
viscosity,  diffiisivity  and  thermal  conduction.  The  models  used  are  discussed  in  detail  in 
[4]  and  [5].  For  the  mixture  viscosity,  p,  the  semi-empirical  formula  of  Wilke  was  used, 
viz.. 


^  =  L  N 


(14) 


where  X;  =  co  /  Mj  is  the  mole  fraction  of  the  i**’  species,  pi  is  the  molecular  viscosity  of 
the  i*  species  at  the  specified  temperature  and  pressure  and 


(b  = 


12 


1  + 


vmJ 


(15) 


1+ 


Ml 

M 


U 


The  individual  species  viscosity,  pi,  is  calculated  from  Chapman-Enskog  kinetic  theory  by 


bii  = 


5 

|7tkM,T 

1671  \ 

1  No 

(16) 
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where  k  is  Boltzmann’s  constant  and  No  is  Avogadro’s  number  and  the  reduced 
temperature,  T*,  is  given  by 


T*  =  — 

Si 


(17) 


where  Oi  and  Sj  are  know  (measured  constants)  for  each  species.  (T*)  is  the  Lennard- 
Jones  potential  for  viscosity. 

The  mixture  thermal  conductivity,  k,  is  calculated  by  a  relationship  similar  to  Eq.  (14) 

(18) 


-N^ 


J=1 


where  Ki  is  the  species  thermal  conductivity  and  is  calculated  from  Euken’s  formula 


K: 


Cpi(T)  + 


5R 
4  Mi 


(19) 


Finally  calculation  of  the  diffusive  fluxes  requires  the  evaluation  of  the  binary  diffusion 
coefficients,  Dy.  In  this  analysis  the  Chapman-Enskog  formulation  is  used 


3  Ttk^N 
8V 


Do  = 


IZ 

r  1 

1  ^ 

JT 

'  H - 

y 

M, 

M- 

V 

v  ^ 

j  / 

po,n„(r) 


(20) 


where 


T’  = 


kT 


(21) 


and  the  Lennard-Jones  parameters,  Oij  and  e,  are  given  by 

Ou=^(o,+o,) 

and 


(22) 


(23) 
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and  Qd  (T*)  is  the  Lennard-Jones  potential  for  diffusion.  It  is  to  be  noted  that  Dy  is 
symmetric,  i.e.,  Dy  =  Dji. 


Multi-component  Diffusive  Flux  Model 

The  diffusive  term  j;  in  Eq.  (4)  consists  of  three  contributions;  (1)  diffusion  due  to 
concentration  gradient  (ordinary  diffosion),  (2)  diffusion  due  to  temperature  gradient 
(Soret  or  thermal  diffusion),  and  (3)  diffusion  due  to  pressure  gradients  (neglected  in  this 
study).  If  only  two  components  are  present,  the  ordinary  diffusive  flux  for  two 
components,  i  and  j,  are  given  by  Pick’s  law. 

j,=-pDyVo.  (24) 

where  Dy  is  the  binary  diffusivity  of  i  into  j.  If  the  carrier  gas  is  the  major  component  and 
all  active  species  are  small  fi-actions,  then  the  mixture  may  be  considered  dilute.  In  this 
case,  each  active  species  may  be  thought  of  as  diffusing  into  the  carrier  gas,  and  the 
corresponding  binaiy  diffusivity  is  an  appropriate  approach.  However,  when  no  dominant 
carrier  gas  is  present  or  when  the  carrier  gas  is  a  light  molecule  such  as  hydrogen  or 
helium,  the  dilute  species  approximation  is  not  appropriate.  It  is  then  necessary  to 
consider  multi-component  mass  transfer.  Multi-component  mass  transfer  is  described  by 
the  Stefan-Maxwell  equations.  Since  solution  of  the  Stefan-Maxwell  equations  is  very 
difficult  and  time  consuming,  in  practice  an  effective  multi-component  diffusivity  for  the  i 
species  into  the  mixture  Di  is  defined.  The  definition  of  D;  is; 


D; 


1-x 


N  Y 

X— 


(25) 


Here  Xi  is  the  mole  fraction.  It  is  noted  that  diflfiisive  fluxes  of  all  species  must  sum  up  to 

N 

zero,  i.e.,  ^  Ji  =  0 ;  only  then  would  all  species  conservation  equations  sum  to  give  the 

i=l 

global  mass  continuity  Eq.  (1).  With  the  form  of  Eq.  (24),  this  condition  is  realizable  in 
practice  if  all  the  diffusion  coefficients  are  equal;  but  this  constraint  is  too  restrictive. 
While  different  approaches  may  be  used  to  satisfy  the  above  constraint,  a  rigorous 
procedure  yields 


M  j=, 


(26) 


as  shown  by  Ramshaw  [6,7].  Though  expression  (26)  is  significantly  more  complex  than 
expression  (24),  it  yields  zero  net  diffusive  flux  for  arbitrary  Di  and  Dj  and  thus  conserves 
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mass.  This  approach  is  better  than  that  often  used  in  the  literature  and  in  those 
commercial  codes  where  the  expression  in  Eq.  (24)  is  used. 

Another  issue  of  concern  when  expression  (24)  is  used  instead  of  (26)  is  that 
usually  one  of  the  species  from  the  set  of  N  is  not  solved,  and  its  mass  or  mole  fraction  is 
obtained  by  subtracting  the  sum  of  the  other  N-1  species  from  unity.  This  is  not  only 
subject  to  round-off  errors,  but  also  may  result  in  an  “asymmetric”  situation,  wherein  the 
final  results  may  depend  on  which  species  is  not  determined  directly  from  its  species 
conservation  equation. 

Thermal  Diffusion  Model 

In  CVD  applications  thermal  (Soret)  diffusion  often  may  be  an  important 
component  of  the  di&sion  process.  The  thermal  diffusion  flux  is  given  by 

ji,T=-D.,TV(lnT)  (27) 

Here,  Dij  is  the  thermal  diffusion  coefficient  for  the  i*  species.  Modelers  generally  rely  on 
empirical  relations  or  data  to  evaluate  Djj.  For  many  species  of  interest  in  structural 
materials  processing,  such  data  are  not  available.  The  present  approach  is  a  first- 
principles-based  procedure  based  upon  the  approach  described  by  Ramshaw  [6,7].  Dij  is 
written  in  the  form 


DiT  =  — M,Ki-o)iY^M:K. 

where  Ki  is  a  thermal  diffusion  coefficient  given  by  the  relationship 


^,=-71;  (A -A,) 

■r  j=\ 


(28) 


(29) 


Here  pij  is  approximated  by  the  expression 


where 


P.i=-P 


4mjYjDy 


JUi- 

2kT 


(30) 


(31) 


and 
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(32) 


'  _2y  "i°# 

'Ci  j=i 

In  these  relations,  k  is  Boltzmann’s  constant,  mj  is  the  mass  of  a  single  molecule  of 
the  i**"  species,  ni  is  the  number  density  of  the  i*  species  and 


m:  +m, 


(33) 


YiYj 

Yi+Yj 


(34) 


(35) 


where  c,  is  the  molecular  collision  diameter  of  the  i*  species  (and  Oy  is  the  total  cross- 
section  for  “i-j”  collisions).  Note  that  t,  represents  the  mean  time  between  collisions  of 
molecules  of  the  i*  species.  The  form  of  Eq.  (30)  differs  from  that  given  in  Ref  12,  and 
was  suggested  by  Ramshaw  [8].  Finally,  note  that  the  form  of  Eq..  (28)  guarantees  that 
the  sum  of  the  thermal  diffusion  fluxes  over  all  species  is  zero. 


Reaction  Models 

In  CVD  applications  two  forms  of  reaction  need  to  be  modeled.  The  first  type  of 
reactions  is  homogeneous  reactions  in  which  the  input  gases  react  to  form  gas  phase 
products.  These  reactions  are  described  with  classical  gas  phase  reaction  models. 
However,  in  CVD  applications  it  is  the  deposition  on  the  wafer  surface  that  provides  the 
desired  end  product.  These  surface  reactions  are  less  well  documented  than  the  gas  phase 
reactions.  Both  types  are  now  discussed. 

Gas  phase  reaction  model 

J 

The  last  term  ^Rij  in  Eq.  (4)  represents  the  sum  of  the  rate  of  production  or 

j=i 

consumption  for  the  i*  species  from  each  reaction).  The  chemical  equation  (conservation 
of  atoms)  for  the  general  elementary  reaction)  can  be  written  as 


i=l 


N 

Zv 

i=l 


X. 


(j  =  1,J) 


(36) 
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where  v|  j  and  v-  j  are  the  stoichiometric  coefficients  of  the  i**"  species  on  the  reactant  and 

product  sides,  respectively,  of  reaction  and  Xi  is  the  chemical  symbol  for  the  i*^  species.  In 
homogeneous  reactions,  the  production  rate  is  given  by 


i.j 


(37) 


where  k,  is  forward  reaction  rate  constant  for  reaction)  which  can  normally  be  expressed 

Ij 

in  terms  of  the  usual  Arrhenius  form,  k,,  is  the  rate  constant  for  the  reverse  reaction,  and 

[Xi]  =  Pj/Mi  is  molar  concentration  of  the  i*^  species.  For  reversible  reactions,  the 
reverse  reaction  rate  constant,  k^,,is  related  to  the  equilibrium  constant  Kc  by 
k,  =  kf  /  .  The  equilibrium  rate  constant  can  be  computed  from  thermodynamic 

properties  (Gibb’s  free  energy)  using  standard  procedures.  Specific  materials  chemistry 
and  reaction  data  to  be  considered  in  this  work  are  discussed  in  a  subsequent  section. 


Surface  process  model 


The  surface  deposition  models  used  in  this  study  are  similar  to  that  described  in 
ref  [9],  The  heterogeneous  reactions  occurring  at  the  interface  between  a  solid  substrate 
surface  and  a  mixture  of  flowing  gases,  the  fluid-substrate  interface,  leading  to  the 
deposition  of  bulk  phase  specie(s)  can  be  viewed  as  a  two  step  process.  In  the  first  step 
the  gaseous  molecules  are  diffused  to  the  surface  where  some  are  adsorbed  onto  the 
various  surface  level  sites.  In  this  study  the  number  of  surface  level  sites  is  assumed  to  be 
a  know  quantity  which  is  a  function  of  the  deposition  surface.  The  second  step  consists  of 
surface  chemical  reactions  and  deposition  of  products  into  the  bulk  phase.  In  general  it  is 
also  possible  to  have  desorption  of  bulk  phase  species  to  the  surface  level  and  of  surface 
species  to  the  gaseous  state.  It  is  clear  from  this  discussion  that  a  surface  site  is  an  active 
participant  in  this  scheme.  In  the  following  discussion  (g)  represents  a  gas  phase  species, 
(s)  denotes  an  adsorbed  surface  species  on  the  top-most  layer  of  the  solid  and  (b)  denotes 
a  deposited  solid.  The  notation  and  handling  procedure  for  surface  species  and  reactions 
can  be  done  in  a  manner  similar  to  the  gas  phase  reactions.  In  analogy  with  Eq.  (36),  the 
surface  reactions  can  be  written 

ivyX,  tvijX,  O-l.D  (38) 

i=l  i=l 

where  in  this  case  K  represents  all  the  species  regardless  of  phase  (gaseous,  surface  and 
bulk),  J*  is  the  total  number  of  surface  reactions  and  X;  is  the  chemical  symbol  of  the  i* 
species.  In  general  all  N  gaseous  species  are  included  in  Eq.  (38).  The  mass  rate  of 
deposition  of  the  i*  species  per  unit  area  per  unit  time,  s, ,  can  be  calculated  from 
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(i=l,K) 


(39) 


J* 

S.  =M,2;Kj-v;,j)qj 

j=l 


where  qj  represents  the  rate  of  progress  variable  for  the  j*  surface  reaction.  In  this  study 
several  forms  of  qj  have  been  used. 


The  most  general  form  for  qj  uses  a  form  similar  to  Eq.  (37) 


^Ij 


(40) 


where  the  same  notation  is  used  as  that  for  Eq.  (37)  except  that  all  variables  are  with 
respect  to  the  surface  reactions.  For  gaseous  species  at  the  surface,  the  concentration  of 
the  i*  species  in  moles  /  cm^  [Xi],  is  defined  the  same  as  for  homogeneous  chemistry.  For 
the  site  species  the  concept  of  a  site  fraction,  Zi,  is  introduced.  The  site  fraction  for  the  i* 
surface  species ,  Z\,  is  related  to  the  concentration  of  the  i*  surface  species,  [Xi],  by 

[X,]  =  ^  (i=N.', . .Nl)  (41) 

CT: 


where  T  is  the  site  density  usually  in  moles  /  cm^  and  ai  is  number  of  sites  occupied  by  the 
i^*  species  and  where  is  the  first  site  species  and  is  the  last  site  species.  Note  that 
the  concentration  of  the  surface  (or  site)  species  is  also  in  moles  /  cm^.  By  definition  the 
sum  of  the  site  fractions  over  all  surface  species  is  equal  to  unity 


n; 


Zz, 

i=N[ 


(42) 


For  bulk  species,  concentrations  are  defined  in  terms  of  the  bulk  species  activities,  ai 

[X,]  =  a,  (i=Nl, . (43) 

where  is  the  first  bulk  species  and  NJ,  is  the  last  bulk  species.  The  activities  are  non- 
dimensional. 

In  this  study  we  have  limited  ourselves  to  cases  of  only  forward  (irreversible) 
reactions  with  only  gaseous  and  site  species  on  the  left  hand  side  of  Eq.  (38).  This 
eliminates  the  need  for  knowledge  of  reverse  reaction  rates  and  bulk  species  activities.  All 
species,  gaseous,  site  and  bulk,  can  occur  on  the  right  hand  side  of  Eq.  (38).  In  the  steady 
state  the  rate  of  production  of  site  species  from  the  gaseous  species  is  equal  to  the  rate  of 
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production  of  the  bulk  species  from  the  site  species.  Therefore,  the  net  production  of 
individual  surface  species  is  zero  or 

s.=0  (i=N^ . ,Ni)  (44) 

Because  there  are  no  bulk  species  on  the  left  hand  side  of  Eq.  (38),  this  results  in  a  system 
of  N5-N5 +1  (in  general)  nonlinear  equations  which  can  easily  be  solved  by  standard 
techniques.  This  method  requires  that  the  system  of  reactions  defined  by  Eq.  (38)  have  all 
site  species  on  both  the  left  and  right  hand  side  of  at  least  one  (not  the  same)  reaction.  If 
this  were  not  the  case  the  only  solution  would  be  a  zero  site  fraction  for  that  species. 
Once  the  system  of  equations  represents  by  Eq.  (44)  is  solved,  the  rate  of  depletion  and 
production  of  the  gaseous  species  and  the  rate  of  production  of  the  bulk  species  can  be 
calculated  by  the  use  of  Eq.  (39).  As  will  be  shown  in  the  Boundary  Conditions  section, 
the  surface  reactions  are  strongly  coupled  to  both  the  fluid  flow  as  described  by  Eqs.  (1)  - 
(4)  and  the  heat  conduction  equation  in  the  solid  surface. 

In  general  the  rate  coefficients,  k^,  and  k,. ,  in  Eq.  (40)  can  be  written  in  the 
Arrhenius  form 


kf  =AiT^^e  RT 


(45) 


where  Aj,  Pj  and  Ej  are  known  coefficients. 

The  sticking  coefficient  model  is  a  subset  of  Eq.  (45)  where  the  relationship  is 


kf,  =v„ 


RT 
271 M, 


(46) 


where  Vi  is  the  “sticking  coefficient”  for  the  i*  species  in  reaction  j.  Sticking  coefficient 
approximations  are  commonly  used  in  surface  chemistry  models.  Although  under  some 
conditions  the  sticking  coefficient  model  may  give  reasonable  results,  it  is  in  general 
unreliable.  For  most  cases  the  surface  chemistry  process  is  more  complex  than  that  which 
can  be  described  by  a  sticking  coefficient  and  hence  has  not  been  used  in  this  study. 

In  some  cases  it  is  desirable  to  use  models  where  k.  and  k.  have  a  non- 
Arrhenius  or  modified  Arrhenius  form  where  the  values  of  k.  and  k  can  be  functions  of 

M  j 

the  concentration  of  the  reacting  species.  In  other  cases  the  values  of  the  qj’s  are 
functions  of  the  concentrations  of  one  or  more  of  the  reacting  species.  In  the  results 
section  the  various  models  will  be  discussed  as  they  are  used. 

Usually  the  desired  end  result  of  a  calculation  is  the  deposition  rate,  f ,  in  cm  /  sec 
of  the  bulk  species  on  the  substrate  surface  is  given  by 
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(47) 


r=y^ 

i^SP. 


where  pi  is  the  density  of  the  i*  bulk  species  and  the  mass  flux  to  the  surface  in  gm  /  cm^  • 
sec,  m ,  is  given  by 


(48) 


Transport  and  Thermodynamic  Properties 

For  the  deposition  of  a  given  material,  the  transport  properties  such  as  viscosity, 
thermal  conductivity  and  difiusivity  must  be  specified  as  a  function  of  temperature  and 
pressure.  These  are  obtained  from  Lennard- Jones  parameters  for  each  of  the  chemical 
species  [10,11],  Note  that  all  physical  properties  vary  across  the  reactor  with  a  variation 
in  temperature  and  composition;  therefore  assumption  of  constant  properties  is  not 
appropriate.  Thermodynamic  properties  such  as  heat  capacity,  enthalpy,  entropy,  and  free 
energy  for  each  species  are  obtained  fi’om  the  NASA  Lewis  database  which  has  been 
added  to  the  code  to  compliment  the  Sandia  National  Laboratory  database.  These 
databases  contains  data  for  many  silane  and  chlorosilane  species  used  in  structural 
materials  preparation. 


Boundary  Conditions 

A  complete  problem  definition  requires  specification  of  boundary  conditions.  At  a 
subsonic  inflow  boundary,  mass  flux  distribution,  flow  angles  and  total  temperature  are 
specified.  This  has  proven  very  effective  for  a  wide  variety  of  cases.  For  species, 
Danckwert's  boundary  condition  is  used  at  the  inlet.  This  condition  equates  the  incoming 
species  mass  fluxes  to  the  net  convective  and  diffusive  fluxes  at  the  inlet.  This  is  an 
alternative  to  specifying  the  mass  or  mole  fi-actions  at  the  inlet.  However,  if  desired  mass 
fractions  may  be  specified.  Outflow  boundary  conditions  require  setting  downstream 
pressure  and  extrapolation  of  species  densities,  temperature  and  velocity. 

At  non-deposition  walls,  the  no-slip  condition  is  imposed  for  velocities.  At  the 
substrate  and  boundaries  where  deposition  occurs,  a  mass  balance  for  each  species  equates 
the  flux  to  the  growth  surface  and  the  net  production  of  that  species  in  surface  reactions. 
In  representing  the  surface  reactions,  one  needs  to  include  detailed  mechanisms  including 
adsorption  and  desorption  kinetics,  as  discussed  earlier  [12].  This  yields 

n(j,+p,U)  =  s,Mi  (49) 
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where  again  the  subscript  i  refers  to  the  i*^  gaseous  species  and  n  refers  to  the  unit  normal 
vector.  The  term  Sj  represents  the  molar  rate  of  deposition  of  species  k  at  the  interface. 
Note  also  that  the  net  growth  rate  results  in  a  non-zero  mass-averaged  normal  velocity  at 
the  growth  surface,  viz., 

N 

n-pU  =  2]s,Mi  (50) 

i=l 

Equation  (50)  can  be  derived  by  summing  Eq.  (49)  over  all  gaseous  species  since  the  sum 
of  all  difRisional  fluxes  must  be  zero.  When  there  is  no  surface  deposition  Eq.  (49)  reverts 
to  the  standard  wall  species  boundary  condition 

nj.=0  (51) 

If  the  effects  of  pressure  and  thermal  diffusion  are  neglected  or  negligible  Eq.  (51) 
becomes  the  simpler  condition 

n-VpOi  =  0  (52) 

In  addition  when  there  is  no  surface  deposition  Eq.  (50)  becomes  the  no  through  flow 
condition 


U  =  0 


(53) 


For  the  thermal  condition  it  is  common  to  specify  either  the  temperature  or  heat 
flux  at  the  reactor  wall  (the  adiabatic  condition  is  a  special  case  of  specified  heat  transfer). 
If  the  heat  flux  is  specified,  it  may  be  necessary  to  include  both  conduction  and  radiation 
effects.  Inclusion  of  radiation  heat  transfer  may  be  necessary  because  of  the  radiation  heat 
exchange  between  heated  susceptor,  walls  and  windows.  For  the  outer  wall,  a  conductive 
description  may  be  adequate.  Although  the  specification  of  the  heat  flux  or  surface 
temperature  is  commonly  used,  it  is  well  known  that  the  substrate  temperatures  are  not 
constant  or  the  heat  transfer  rates  are  imprecisely  known.  A  more  general  approach  is  to 
solve  the  heat  conduction  equation  in  the  solid  simultaneously  with  the  fluid  dynamic 
equations.  In  this  case  it  is  necessary  to  specify  a  compatibility  condition  at  the  between 
the  fluid  and  the  solid  substrate. 


When  the  heat  conduction  equation  is  solved  in  conjunction  with  the  fluid 
mechanical  equations  (  the  so-called  conjugate  technique),  the  energy  transfer  within  the 
solid  material  is  governed  by  the  heat  conduction  equation  of  the  form 


Ps^D 


dt 


V  k  VT +Q, 


(54) 
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where  ps  k,  and  T,  are  the  density,  specific  heat,  thermal  conductivity  and  temperature 

respectively  of  the  solid.  Qs  represents  any  non-conductive  heat  generation  within  the 
solid  as  for  example  might  be  due  to  inductive  heating.  The  solid  material  properties  do 
not  have  to  be  constant  and  are  often  represented  as  a  functions  of  temperature.  The 
compatibility  condition  at  the  interface  of  the  fluid  and  the  solid  surface  ensures  that  the 
temperature  and  the  heat  transfer  across  the  interface  are  consistent.  Therefore,  at  the 
interface  junction,  the  temperature  of  the  fluid  and  the  solid  must  be  the  same 

Tg  =  X 

while  the  energy  balance  at  the  interface  is  given  by 

"•S'^Tg-ZWi(T)  =  QR+n-K,VT,  (56) 

i=l 

The  terms  on  the  left  hand  side  of  Eq.  (56)  represent  the  conductive  heat  transfer  from  the 
fluid  and  the  heat  transfer  due  to  deposition.  Species  are  summed  over  all  K  (gaseous, 
solid  and  bulk)  species  at  the  interface.  The  terms  on  the  right  hand  side  of  the  equation 
represent  radiation,  Qr,  and  conductive  heat  transfer  into  the  solid.  In  this  equation  the 
subscript  g  represents  the  gaseous  (fluids)  values  at  the  interface.  In  practice  the  fluid 
dynamics  equations  and  the  heat  conduction  equation  are  solved  independently  with  the 
above  compatibility  equations  enforced  as  the  boundary  conditions. 

When  the  conjugate  problem  is  considered  in  this  study,  the  solution  procedure  is 
to  first  solve  the  fluid  dynamic  equations  with  specified  interface  temperatures.  The  heat 
transfer  rate  from  the  fluids  to  the  solid  surface  and  the  deposition  energy  transfer  (the 
first  and  second  terms  of  the  left  hand  side  of  Eq.  (56)  are  calculated  and  then  imposed  on 
the  heat  conduction  equation  as  the  boundary  condition.  The  heat  conduction  is  then 
solved,  the  surface  temperatures  calculated  and  imposed  on  the  fluid  dynamics.  The 
overall  process  is  repeated  until  both  the  fluid  dynamic  equations  and  the  heat  conduction 
equation  are  converged. 

Solution  and  Numerical  Procedure 

General 


The  basic  numerical  procedure  and  in-house  computational  fluid  dynamics  code 
has  been  developed  by  SRA  personnel  over  many  years.  The  current  capability  includes 
two  and  three  dimensional  analysis,  reacting  flow,  two  phase  flow  and  multiple  turbulence 
models.  More  recently,  as  will  be  discussed  subsequently,  SRA  under  its  own  funds  has 
developed  an  incompressible  (or  very  low  Mach  number)  capability  in  which  pressure 
coefficient  appears  specifically  as  a  dependent  variable.  Applications  of  the  code, 
numerical  procedure  and  grid  capability  have  been  made  to  a  wide  range  of  practical 
problems  for  DOD  and  NASA  including  turbo-machinery,  primary  and  secondary  gas 
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paths,  re-entry  vehicles,  rocket  motors,  inlets,  seals,  chemical  lasers,  combustion,  two  and 
multi-phase  flows,  crystal  growth  and  CVD.  Examples  can  be  found  in  [12-23], 


Algorithm 

The  numerical  algorithm  which  forms  the  basis  for  this  code  is  described  below. 
The  procedure  to  solve  the  governing  equations  is  a  consistently  split  linearized  block 
implicit  scheme  originally  developed  by  Briley  and  McDonald  [24]  at  SRA  and  embodied 
in  a  computer  code  termed  MINT.  The  basic  algorithm  has  been  further  developed  and 
applied  to  both  laminar  and  turbulent  flows  (see  Briley  and  McDonald  [25]).  The  method 
can  be  outlined  as  follows:  the  governing  equations  in  a  non-dimensional  form  are 
replaced  by  an  implicit  time  difference  approximation.  Terms  involving  nonlinearities  at 
the  unknown  time  level  are  linearized  by  Taylor  series  expansion  about  the  solution  at  the 
previous  known  time  level,  and  spatial  difference  approximations  are  introduced.  The 
result  is  a  system  of  multi-dimensional  coupled  (but  linear)  difference  equations  for  the 
dependent  variables  at  the  unknown  or  implicit  time  level.  To  solve  these  difference 
equations,  the  Douglas-Gunn  procedure  for  generating  alternating  direction  implicit  (ADI) 
splitting  schemes  is  used.  This  ADI  splitting  technique  leads  to  systems  of  coupled  linear 
difference  equations  having  narrow  block-banded  matrix  structures  which  can  be  solved 
efficiently  by  standard  block-elimination  methods. 

Low  Mach  number  or  incompressible  flow  considerations 

Very  low  Mach  number  or  incompressible  flow  can  lead  to  stability  problems  for 
the  solution  procedures  described  above.  Basically  at  very  low  Mach  numbers  the 
solution  matrix  becomes  ill-conditioned  leading  to  slow  or  non-convergence  and,  perhaps, 
to  loss  of  accuracy.  The  difficulties  increase  when  significant  temperature  variation  is 
present  as  is  normally  the  case  in  CVD  processes. 

Under  a  corporate  IR&D  program,  SRA  has  extended  its  Navier  Stokes  solution 
procedure  to  low  Mach  number  or  incompressible  flows.  With  this  approach  density,  p,  is 
replaced  by  the  pressure  coefficient,  Cp,  as  a  dependent  variable.  The  global  continuity, 
momenta  and  species  continuity  governing  partial  differential  equations,  Eqs,  (1)  -  (4),  are 
recast  in  the  form 


a 


+  V-pU  =  0 


(57) 


^  +  VpUU  =  -ip„UiVc,+VT+pg  (58) 

^  +  VphU  =  -V.q  +  ip.U'^+4>  +  Q„  (59) 

Ot  2  Dt 
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and 


d\ 
a  — 


pQ^jCpj 

“St 


+  v-p0,u  =  -v.ji+XRij 

j=i 


(60) 


where  a  is  an  arbitrary  constant  and  the  pressure  coefficient,  Cp,  is  defined  by 


Cp  = 


P-P. 

ip.U’. 


(61) 


For  solutions  where  only  steady  state  is  of  interest  the  solution  obtained  from  Eqs. 
(57)  -  (60)  are  the  same  as  the  solution  obtained  fi-om  Eqs.  (1)  -  (4).  The  main  advantage 
in  solving  this  modified  set  of  equations  is  that  the  dependent  variable,  Cp,  does  not  have 
the  large  spatial  variation  that  the  density,  p,  has  especially  when  there  are  large 
temperature  gradients.  Hence  the  convergence  and  stability  characteristics  of  the 
governing  equations  are  better  behaved.  When  Eq.  (5)  is  substituted  into  Eq.  (61)  a  new 
form  of  the  equation  of  state  is  obtained 


P  = 


M 

RT 


P„  + 


(62) 


Because  of  the  artificial  time  terms  in  the  global  continuity  and  the  species 
conservation  equations  the  above  equations  cannot  be  used  for  time  accurate  solutions. 
Time  accurate  forms  of  Eqs.  (1)  -  (4)  can  be  obtained  by  augmenting  the  physical  time 

derivatives  with  artificial  time  terms  (as  described  above)  of  the  form  —  where  x 

ox 

represents  an  artificial  inner  iteration  time  term  chosen  to  accelerate  inner  convergence. 
For  this  case  the  governing  equations  take  the  form 


^  P£p 

^+a  ^  ^  ■i  +  V-pU  =  0 
at  dx 


■^  +  ^+V-pUU  =  4p«UiVCp+V-x+pg 
dt  dx  2 


M  +  M  +  V-phU^-V-q  +  lp^Ui^  +  O  +  Q, 
at  ax  2  Dt 


(63) 

(64) 

(65) 
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and 


at 


+  a 


/ 

d 

V 


2  J 


5t 


+vp®,u=-v.j,+2;R, 

j=i 


(66) 


Eqs.  (63)  -  (66)  are  the  same  as  Eqs.  (1)  -  (4)  if  the  artificial  time  terms  are  neglected. 
Time  accurate  solution  of  Eqs.  (63)  -  (66)  are  obtained  by  using  an  inner  iteration  -  outer 
iteration  solution  procedure.  In  the  inner  iteration  procedure  the  time  step  t  is  chosen  to 
optimize  inner  convergence  of  the  iteration  procedure.  The  time  accurate  portion  of  the 
solution  is  determined  fi-om  the  physical  time  term. 
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•RESULTS 


General 

The  analysis  described  above  has  been  applied  to  a  number  of  test  cases.  Although 
the  physics,  homogeneous  and  heterogeneous  chemistries  and  geometries  of  the  cases  are 
significantly  different,  there  are  certain  common  features  that  apply  to  the  general  form  of 
the  input  and  output  of  all  the  cases.  The  cases  have  all  been  run  with  the  MINT 
computer  code  which  solves  the  Navier-Stokes  equations  in  the  gaseous  flow  and  the  heat 
conduction  equation  (if  desired)  in  the  solid  substrate  with  the  appropriate  compatibility 
conditions  and  heterogeneous  surface  chemistry  applied  at  the  interface.  The  input  for  all 
cases  is  outlined  in  Table  I.  Input  is  divided  into  three  parts:  (1)  Reactor  Input,  (2) 
Substrate  Input  and  (3)  Deposition  Surface  Input.  The  reactor  and  substrate  geometries 
must  be  compatible,  i.e.,  at  the  fluid-substrate  interface  the  two  geometries  must  be 
coincident.  The  grids  do  not  have  to  match  grid  point  to  grid  point.  Interpolation 
procedures  can  be  used  to  transfer  information  from  the  fluids  grid  to  the  substrate  grid 
and  vice  versa.  For  all  cases  considered,  the  flow  is  subsonic  with  the  Mach  number  being 
on  the  order  of  10'^  to  10"^.  This  requires  the  use  of  the  subsonic  inlet  and  exit  boundary 
conditions  described  in  the  boundary  conditions  section  and  outlined  in  Table  I.  The 
homogeneous  chemistry  used  will  be  described  in  detail  for  each  case.  In  some  cases  the 
homogeneous  chemistry  involves  many  chemical  species  with  numerous  chemical  reaction 
described  by  Arrhenius  type  reactions.  In  still  other  cases  the  homogeneous  chemistry 
uses  a  “global”  description  where  a  semi-empirical  formula  is  used  to  describe  the 
homogeneous  chemistry.  In  other  cases  there  are  no  homogeneous  reactions.  All  the 
reactions  take  place  at  the  fluids-substrate  interface  in  the  form  of  heterogeneous 
reactions. 

Substrate  conditions  are  usually  straight-forward.  In  the  temperature  range  of 
interest  the  density,  specific  heat  and  thermal  conductivity  of  the  substrate  material  are 
normally  considered  to  be  known  constants.  Usually  in  CVD  applications  there  is  a 
desired  target  temperature  at  the  fluids-substrate  interface.  In  practice  a  thermocouple  is 
embedded  in  the  surface  to  approximately  monitor  this  temperature  and  the  amount  of 
energy  applied  to  the  substrate  is  modified  until  that  desired  temperature  is  achieved.  In 
the  numerical  simulation  the  surface  temperature  is  either  directly  set  (if  the  conjugate 
technique  is  not  used)  or  the  surface  temperature  is  numerically  monitored  (if  the 
conjugate  technique  is  used)  and  the  substrate  backface  boundary  condition  is  modified 
until  the  desired  temperature  is  achieved. 

The  input  for  the  heterogeneous  chemistry  is  similar  to  the  homogeneous 
chemistry.  Both  Arrhenius  and  non-Arrhenius  heterogeneous  chemistries  have  been  used 
as  well  as  “global”  techniques.  When  the  Arrhenius  site  fraction  techniques  have  been 
used  the  site  or  surface  species  deposition  rates  are  determined  from  the  equilibrium 
condition  previously  described.  Therefore  it  is  not  necessary  to  know  the  surface  species 
thermo-chemical  properties  necessary  to  calculate  the  species  enthalpy  in  Eq.  (56) 
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REACTOR  INPUT 


Reactor  Geometry 
Boundary  Conditions 

Input  Flow  Conditions 
Mass  Flux 

Stagnation  Temperature 
Composition 
Exit  Conditions 

Specified  Static  pressure 
Symmetry  Conditions 
Wall  or  Surface  Conditions 
Deposition 
Heat  Transfer 

Homogeneous  Rate  Chemistry 
Gaseous  Species  Thermo-chemical  Data 

SUBSTRATE  INPUT 

Substrate  Geometry 
Material  Properties 
Surface  Boundary  Conditions 
Deposition 
Symmetry 

Heat  Transfer  Specified 
Temperature  Specified 

DEPOSITION  SURFACE  INPUT 

Heterogeneous  Rate  Chemistry 

Surface  and  Bulk  Species  Thermo-chemical  Data 

Site  Density 


Table  I.  CVD  Reactor  Input 

The  form  of  the  output  obtained  for  all  cases  is  shown  in  Table  II.  In  the  reactor 
all  flow  variables,  dependent  as  well  as  derived  variables,  are  produced  as  part  of  the 
output  and  can  be  displayed,  plotted  or  averaged  as  desired  by  the  user.  Species 
distribution  can  be  displayed  as  mass  fractions,  mole  fractions  or  concentrations.  Flow 
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REACTOR  OUTPUT 


Flow  Patterns 
Temperature  Field 
Pressure  Field 
Species  Distribution 
Other  Variables 

SUBSTRATE  OUTPUT 

Temperature  Field 

DEPOSITION  SURFACE  OUTPUT 

Film  Growth  Rate 

Removal  of  Gases 
Deposition  of  Bulk  Species 
Growth  Rate  Uniformity 
Film  Composition 


Table  II.  CVD  Reactor  Output 

patterns  can  be  obtained  to  aid  in  the  design  of  more  efficient  CVD  systems.  When  the 
conjugate  technique  is  used,  the  temperature  distribution  and  heat  transfer  rates  in  the 
substrate  can  be  displayed.  The  ultimate  objective  of  the  CVD  analysis  is  to  obtain 
information  about  the  film  deposition  process.  One  can  display  film  growth  rates, 
uniformity  and  composition  as  a  function  of  the  position  on  the  fiuids-substrate  interface. 

Test  Cases 

A  series  of  six  test  cases  have  been  run  in  this  study.  These  cases  have  been 
chosen  for  a  variety  of  reasons.  Some  were  chosen  because  of  the  need  to  obtain 
information  about  critical  DOD  CVD  programs,  while  others  were  chosen  to  demonstrate 
the  ability  of  the  above-described  technique  to  duplicate  more  academic  experimental  data. 
All  cases  were  run  to  the  steady  state  using  Eqs.  (57)  -  (60)  as  the  governing  partial 
differential  equations.  Two  of  the  six  case  were  run  with  the  conjugate  technique, 
duplicating  the  calculations  performed  for  the  same  physical  device  when  using  the 
specified  temperature  interface  condition,  the  non-conjugate  mode.  The  presentation  of 
each  case  consists  of  information  presented  in  a  similar  format.  This  includes;  (1)  a 
schematic  drawing  and  brief  verbal  description  of  the  physical  device,  (2)  a  table  that 
contains  flow  input  information  necessary  to  perform  the  numerical  calculation  (3)  tables 
showing  the  homogeneous  and  heterogeneous  chemical  reactions  and  associated  rate 
constants  and  (4)  a  verbal  and  graphical  description  of  the  results  and  what  the  results 
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demonstrate.  When  a  large  number  of  cases  were  run,  as  was  the  case  for  the  first  four 
test  cases,  this  section  is  described  in  the  written  text. 

Silane-Propane  and  Methyltrichlorosilane  -  Cases  I  -  IV 

The  first  four  cases  considered  in  this  study  use  the  same  basic  CVD  reactor 
geometry.  A  schematic  of  the  device  is  shown  in  Fig.  1.  It  represents  a  typical  vertical 
CVD  reactor  employed  in  semiconductor  processing. 


! 


I 

I 


Fig.  1 .  Schematic  Diagram  for  CVD  Reactor 

Flow  containing  the  precursor  gas  and  the  carrier  gas  for  SiC  deposition  enters 
from  a  surrounding  inlet  section  on  the  top  of  the  device,  proceeds  towards  the  centerline 
and  is  ducted  down  a  metal  pipe  onto  the  heated  susceptor  where  the  deposition  occurs. 
Flow  then  proceeds  in  a  circuitous  manner  away  from  the  centerline  eventually  exiting  at 
the  bottom  of  the  device  as  shown.  The  reactor  geometry  is  cylindrically  symmetric  and 
hence  the  computations  are  two-dimensional  axisymmetric,  i.e.,  the  swirl  component  of 
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velocity  is  zero  and  there  is  no  angular  variation  of  any  of  the  flow  variables.  Critical 
dimensions  are  noted  on  Fig.  1 .  dp  represents  the  diameter  of  the  pipe,  Zp  represents  the 
distance  from  the  susceptor  surface  to  the  end  of  the  pipe  and  d,  is  the  deposition  surface 
diameter.  Initially  values  of  dp  and  Zp  were  varied  to  determine  ‘optimum’  values  for 
deposition  rate  and  consistency  of  thickness.  These  values  were  found  to  be  2  cm  and  1 
cm  respectively.  The  value  of  ds  was  always  take  as  5  cm.  The  temperature  of  the 
susceptor  (substrate)  was  varied  either  directly  (by  specifying  the  susceptor  temperature) 
or  indirectly  (by  solution  of  the  heat  conduction  equation  in  the  susceptor  and  appropriate 
back  face  temperature)  and  the  outer  walls  were  assumed  to  be  at  room  temperature.  The 
inlet  gas  was  at  290  K.  For  these  cases,  two  feed  precursor  gases  were  considered;  (1) 
Silane-Propane  and  (2)  Methyltrichlorosilane  (MTS).  Both  of  the  feed  gases  used 
hydrogen  as  a  carrier  gas.  The  molar  ratio  of  the  precursor  gas  to  the  hydrogen  carrier 
gas  considered  were  10%,  20%  and  40%  with  the  majority  of  calculations  performed  at 
10%.  .  The  homogenous  and  heterogeneous  reactions  sets  used  in  these  calculations  are 
shown  in  Tables  III  &  IV  for  Silane-Propane  and  Tables  V  &  VI  for  MTS.  The 
homogenous  Silane-Propane  chemistry  considered  13  individual  species  while  the  MTS 
homogenous  chemistry  considered  15  species.  The  heterogeneous  Silane-Propane  surface 
reactions  require  8  additional  chemical  species  while  the  corresponding  MTS  surface 
reactions  require  7  additional  species.  Thus  the  Silane-Propane  cases  required  the  solution 
of  17  coupled  partial  differential  equations  and  the  MTS  case  required  the  solution  of  19 
coupled  partial  differential  equations.  The  finite  difference  grid  used  for  the  computations 
consisted  of  58  grid  points  in  the  axial  direction  and  45  grid  points  in  the  radial  direction. 
Critical  regions,  such  as  the  substrate  and  wall  regions,  used  a  concentration  of  grid  cells 
in  these  regions  relative  to  other  places,  the  basic  philosophy  being  that  a  dense 
concentration  of  cells  is  used  where  gradients  are  large. 

When  one  considers  the  design  of  a  CVD  device  many  factors  must  be  considered. 
First  goodness  criteria  must  be  established  by  which  you  can  evaluate  the  desired  physical 
results.  For  this  case  two  arbitrary  criteria  were  used.  The  first  is  the  deposition  rate 
(usually  measured  in  nm/min).  The  second  is  the  uniformity  of  the  deposition  rate.  Large 
growth  rates  are  considered  to  be  desirable  while  a  large  radial  variation  was  considered  to 
be  undesirable  and  a  more  uniform  variation  was  considered  to  be  more  desirable.  Design 
parameters  that  one  would  consider  certainly  would  include  at  least  the  following:  (1)  the 
mixture  ratio  of  delivery  precursor  gas  to  carrier  gas,  (2)  the  total  flow  rate  of  precursor 
and  carrier  gas,  (3)  the  delivery  pressure  of  the  gas  and  (4)  the  temperature  of  the 
susceptor  surface.  Also  one  could  consider  the  above-mentioned  dimensions  of  the  CVD 
device.  If  all  these  possibilities  were  taken  into  consideration  the  matrix  of  possibilities 
would  require  a  very  large  number  of  runs;  therefore  some  compromise  needed  to  be 
made.  The  technique  used  was  to  consider  the  geometry  as  frozen  and  to  choose  flow 
rate,  pressure  and  susceptor  pressure  at  set  values  (it  will  later  be  shown  that  these  were 
not  exactly  arbitrary  values).  The  values  chosen  for  total  volumetric  flux  of  precursor  and 
carrier  gas  was  set  as  600  standard  cubic  centimeters  per  minute  (seem).  The  notation 
seem  refers  to  the  amount  of  volumetric  flow  rate  that  would  correspond  to  the  mass  flux 
that  occur  at  reference  conditions  of  1  atmosphere  and  room  temperature.  The  pressure 
was  chosen  as  15.2  torr  and  the  temperature  of  the  susceptor  was  chosen  as  1300  K.  The 
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Reactions 

_ _ 

Ej 

1) 

2H 

+ 

M 

0 

H 

+ 

M 

2.95x10** 

-1.0 

0 

2) 

CH, 

+ 

H 

<=> 

CH 

4 

H 

1.259x10*^ 

0.0 

11900 

3) 

CH4 

+ 

M 

<=> 

CH 

+ 

H 

4- 

M 

1.413x10*'* 

0.0 

88400 

4) 

2CH3 

<=> 

C2H 

8.913x10*^ 

0.0 

0 

5) 

C2H6 

+ 

H 

0 

C2H 

4- 

H 

5.40x10^ 

3.5 

5210 

6) 

C2H5 

4- 

M 

0 

C2H 

+ 

H 

+ 

M 

1.99x10*^ 

0.0 

30000 

7) 

2CH3 

0 

C2H 

H 

2.80x10*^ 

0.0 

13592 

8) 

C2H6 

-(- 

CH 

<=> 

C2H 

CH 

5.50x10* 

4.0 

8300 

9) 

C2H4 

-h 

H 

<=> 

C2H 

2.21x10*^ 

0.0 

2066 

10) 

C2H4 

+ 

M 

0 

C2H2 

+ 

H 

M 

1.50x10*^ 

0.0 

55800 

11) 

C3Hg 

0 

C2H 

+ 

CH 

1.70x10**^ 

0.0 

84840 

12) 

SiH, 

0 

SiH2 

4- 

H 

6.67x10^^ 

-4.795 

63450 

13) 

SiH 

0 

SiH3 

+ 

H 

3.69x10** 

0.0 

93000 

14) 

SiH 

-h 

H 

<=> 

SiH3 

+ 

H 

1.46x10** 

0.0 

2500 

15) 

SiH2 

H 

0 

SiH 

4- 

H 

1.39x10** 

0.0 

2000 

16) 

SiH2 

+ 

H 

SiH3 

3.81x10** 

0.0 

2000 

17) 

SiH 

4- 

H 

SiH3 

3.45x10** 

0.0 

2000 

Table  III.  Homogeneous  Reactions:  SiHt  -  CsHg  -  H2  System 


Reactions 


A 


n. 


1) 

2) 

3) 

4) 

5) 

6) 

7) 

8) 

9) 

10)  SiH 

11) 

12) _ 


CH3 

CH4 

C2H5 

C2H4 

C2H2 

SiH2 

SiH4 

SiHs 


+  Si(s) 

+  Si(s) 

+  2Si(s) 

+  2Si(s) 

+  2Si(s) 
2CH(s) 
+  C(s) 

+  C(s) 

+  C(s) 

+  C(s) 
2SiH(s) 
SiH2(s) 


H2 

2H2 

2H2 

2H2 

H2 

H2 

H2 

H2 

H2 

H2 


CH(s) 

C(s) 

C(s) 

2C(s) 

2C(s) 

2C(s) 

SiH2(s) 

SiH2(s) 

SiH(s) 

SiH(s) 

2Si  (s) 

Si  (s) 


+  Si(b) 

+  Si(b) 

+  CH(s)  +  2Si(b) 
+  2Si(b) 

+  2Si(b) 

+  C(b) 

+  C(b) 

+  C(b) 

+  C(b) 


8.67x10" 
4.20x10’ 
5.76x10^® 
9.37x10" 
1.22x10'^ 
2.25x10’'' 
6.12x10" 
3.18x10'® 
6.03x10" 
6.23x10" 
2.25x10’'' 
2.91x10 


14 


0.5 

0.5 

0.5 

0.5 

0.5 

0.0 

0.5 

0.5 

0.5 

0.5 

0.0 

0.0 


0 

0 

0 

0 

0 

61000 

0 

18678 

0 

0 

61000 

9000 


Table  IV.  Heterogeneous  Reactions:  SiHt  -  CsHg  -  H2  System 

first  variable  to  be  investigated  was  the  molar  ratio  of  precursor  gas  to  carrier  gas  for 
Silane-Propane.  Three  molar  ratios  of  10%,  20%  and  40%  were  considered.  Fig.  2 
shows  the  results.  As  can  be  seen  in  Fig.  2,  the  radius  varies  fi-om  -2.5  cm  to  2.5  cm  or  a 
total  distance  of  ds  =  5.0  cm.  The  deposition  rate  is  inversely  proportional  to  the  molar 
ratio.  Also  the  uniformity  of  deposition  remains  relatively  constant  with  respect  to  radius 
with  some  slight  variation  at  a  ratio  of  10%.  Therefore  using  our  above  criteria  of 
goodness,  it  was  decided  to  use  a  molar  ratio  of  10%  for  all  subsequent  investigations.  It 
is  recognized  that  this  decision  has  a  certain  amount  of  arbitrariness  and  one  might  even 
make  an  argument  for  a  molar  ratio  of  0.2  based  on  a  slightly  more  uniform  deposition. 
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Reactions 

Ai 

Pi 

Ej 

1) 

2H 

+  M 

0  H2 

+ 

M 

2.95x10** 

-1.0 

0 

2) 

CH4 

+  H 

0  CH3 

+ 

H2 

1.259x10*'* 

0.0 

11900 

3) 

CH, 

+  M 

0  CH3 

+ 

H 

+ 

M 

1.413x10*’ 

0.0 

88400 

4) 

2CH3 

0  C2H6 

8.913x10*^ 

0.0 

0 

5) 

C2H6 

+  H 

C2H5 

+ 

H2 

5.40x10^ 

3.5 

5210 

6) 

C2H5 

+  M 

0  C2H4 

+ 

H 

M 

1.99x10*^ 

0.0 

30000 

7) 

2CH3 

0  C2H5 

H 

2.80x10*^ 

0.0 

13592 

8) 

C2H6 

+  CH3 

0  C2H5 

+ 

CH4 

5.50x10* 

4.0 

8300 

9) 

C2H4 

+  H 

0  C2H5 

2.21x10*^ 

0.0 

2066 

10) 

C2H4 

+  M 

0  C2H2 

+ 

H2 

+ 

M 

1.50x10*^ 

0.0 

55800 

11) 

Si(CH3)Cl3 

•0  SiCb 

CH3 

7.60x10*'* 

0.0 

69312 

12) 

SiHCl3 

+  H 

<»  SiCb 

+ 

H2 

2.47x10*^ 

0.0 

2543 

13) 

SiHClsH 

0  SiCl2 

+ 

HCl 

2.60x10** 

0.0 

47000 

14) 

SiCU 

+  H 

0  SiCb 

+ 

HCl 

1.50x10*^ 

0.0 

3400 

Table  V.  Homogeneous  Reactions:  MTS  -  H2  System 


Reactions 


tL. 


1) 

2) 

3) 

4) 

5) 

6) 

7) 

8) 

9) 

10)  H 


CHj 

CH4 

C2H5 

C2H4 

C2H2 

SiCb 

H2 

SiCb 


+  Si(s)  =>  H2  +  CH(s)  +  Si(b) 

+  Si(s)  ^  2H2  +  C(s)  +  Si(b) 

+  2Si(s)  =>  2H2  +  C(s)  +  CH(s) 

+  2Si(s)  =>  2H2  +  2C(s)  +  2Si(b) 

+  2Si(s)  =>  H2  +  2C(s)  +  2Si(b) 

2CH(s)  =>  H2  +  2C(s) 

+  C(s)  =j>  SiCbCs)  +  C(b) 

+  SiCbCs)  =>  2HC1  +  Si{s) 

+  C(s)  =>  SiCl3(s)  +  C(b) 

+  SiCbts)  =>  HCl  +  SiCbCs) _ 


+  2Si(b) 


8.67x10" 

4.20x10’ 

5.76x10^° 

9.37x10" 

1.22x10" 

2.25x10^" 

3.38x10" 

6.22x10'* 

2.90x10" 

6.22x10’ 


0.5 

0.5 

0.5 

0.5 

0.5' 

0.0 

0.5 

1.0 

0.5 

1.0 


0 

0 

0 

0 

0 

61000 

0 

0 

0 

0 


Table  VI.  Heterogeneous  Reactions:  MTS  -  H2  System 


The  second  variable  considered  was  the  overall  flow  rate  for  both  Silane-propane  and 
MTS.  The  results  are  shown  in  Figs.  3-4  respectively.  In  both  cases  it  can  clearly  be  seen 
that  600  seem  is  the  best  flow  rate.  It  should  also  be  noted  that  if  a  different  geometric 
configuration  had  been  chosen  the  results  might  have  been  different  (not  having  performed 
those  calculations  this  cannot  be  said  with  a  complete  degree  of  certainty). 

The  third  variable  considered  was  the  pressure.  Three  pressures,  15.2  torr,  38.0 
torr  and  76.0  torr  were  chosen  as  representative  of  CVD  applications.  Figs.  5-6  show  the 
results  for  injector  pressures  of  Silane-Propane  and  MTS  respectively.  It  is  clearly  seen 
that  at  76.0  torr  there  is  a  large  radial  variation  of  growth  rate.  For  both  cases  at  15.2  torr 
the  growth  rate  is  relatively  uniform.  As  can  be  seen  at  38.0  torr  the  variation  in  growth 
rate  is  only  slightly  less  uniform  than  at  15.2  torr.  However,  a  pressure  of  15.2  torr  was 
chosen  as  the  reference  condition  for  the  remainder  of  this  study. 
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Fig  2.  Silane-Propane  Growth  Rate  vs.  Radial  Distance  as  a  Function  of  Molar 

Ratio 


Fig  3.  Silane-Propane  Growth  Rate  vs.  Radial  Distance  as  a  Function  of 
Volumetric  Flow  Rate 
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Fig  4.  MTS  Growth  Rate  vs.  Radial  Distance  as  a  Function  of  Volumetric  Flow 

Rate 
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Fig.  5.  Silane-Propane  Growth  Rate  vs.  Radial  Distance  as  a  Function  of  Inlet 

Pressure 
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Fig.  6.  MTS  Growth  Rate  vs.  Radial  Distance  as  a  Function  of  Inlet  Pressure 
The  fourth  variable  considered  was  the  susceptor  temperature  shown  in  Figs.  7-8. 


Fig.  7.  Silane-Propane  Growth  Rate  vs.  Radial  Distance  as  a  Function  of  Susceptor 

Temperature 
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Three  temperatures  were  chosen  for  the  Silane-Propane  case,  1300  K,  1450  K  and  1600  K 
and  three  separate  temperatures  were  chosen  for  the  MTS  case,  1000  K,  1150  K  and  1300 
K.  These  temperatures  were  chosen  to  be  representative  of  CVD  applications.  For  the 
Silane-Propane  case,  the  1000  K  case  yielded  a  value  too  small  to  appear  on  the  figure. 
Even  the  1 150  K  case  yields  a  small  deposition  rate.  The  1300  K  case  yields  a  significantly 
larger  deposition  rate  with  a  fairly  uniform  radial  variation.  For  the  MTS  case  it  is  clear 
that  the  1300  K  case  yields  the  most  uniform  deposition.  At  1450  K  and  1600  K  the  radial 
variation  has  no  region  of  constant  deposition  rate. 


Fig.  8.  MTS  Growth  Rate  vs.  Radial  Distance  as  a  Function  of  Susceptor 

Temperature 

To  test  the  conjugate  capability  two  test  cases  were  run,  one  with  a  Silane- 
Propane  precursor  and  one  with  an  MTS  precursor.  Both  cases  were  run  at  the  above 
derived  ‘optimum’  or  reference  conditions  of  600  seem  of  precursor  and  carrier  gas,  10% 
precursor  gas,  at  an  inlet  pressure  of  15.2  torr  and  a  nominal  susceptor  temperature  of 
1300  K.  The  susceptor  temperature  was  set  indirectly  by  setting  the  backface  temperature 
of  the  substrate  at  a  nominal  temperature.  A  converged  solution  was  obtained  and  the 
average  susceptor  temperature  calculated.  Modifications  were  then  made  to  the  backface 
temperature  until  the  desired  or  target  susceptor  temperature  was  obtained.  This  can  be 
viewed  as  a  simulation  of  a  resistance-heating  coil  placed  on  the  back  of  the  substrate 
pedestal.  One  could  then  vary  the  current  until  the  desired  susceptor  temperature  was 
obtained.  The  geometry  is  the  same  as  shown  in  Fig.  1 .  The  substrate  was  assumed  to  be 
made  from  stainless  steel  with  a  specific  gravity  of  7.9,  a  specific  heat  of  594.4  joules  /  kg- 
K,  and  a  thermal  conductivity  of  16.3  joules  /  m-sec-K.  For  the  Silane-Propane  case  the 
backface  of  the  substrate  was  set  to  2378  K  and  for  the  MTS  case  the  backface  of  the 
substrate  was  set  to  2320  K.  In  both  cases  this  led  to  a  temperature  of  approximately 
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1300  K  on  the  susceptor/fluids  interface.  To  allow  for  a  smooth  variation  of  temperature 
at  the  lower  left  and  right  sides  of  the  substrate,  the  temperature  in  these  regions  was 
smoothly  lowered  to  290  K  over  the  last  0.84  cm  of  the  substrate.  Fig.  9  shows  the 
temperature  contours  for  both  cases.  This  figure  shows  only  the  right-hand-side  of  the 
device  portrayed  in  Fig.  1.  The  substrate  is  outlined  with  a  black  box  in  the  lower  left- 
hand-side  of  both  cases.  The  0.84  cm  region  where  the  backface  temperature  is  lowered 
to  290  K  is  evident  in  this  figure.  It  is  important  to  note  the  continuity  of  temperature 
across  all  substrate/fluid  boundaries,  the  so-call  compatibility  condition.  Not  oriy  the 
temperature  but  also  the  heat  transfer  is  continuous  across  these  boundaries.  While  it  is 
hard  to  discern  fi-om  the  figure,  the  heat  transfer  across  the  MTS/substrate  boundary  is 
slightly  larger  than  for  the  comparable  SUane-Propane  case.  This  is  probably  due  to  the 
difference  in  the  thermal  conditions  of  the  precursor  gases.  As  might  be  expected  the 
deposition  rates  were  virtually  the  same  as  for  the  non-conjugate  runs 


Silane  Propane  MTS 


0.0  1100.0  2200.0  OjO  1100.0  2200.0 


Figure  9.  Temperature  Contours  for  Silane-Propane  and  MTS  at  Standard  Conditions 
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The  purpose  of  the  above  exercise  -was  to  demonstrate  the  ability  of  the  analysis  to 
help  the  engineer  to  design  a  CVD  device  that  can  in  some  sense  be  optimized  for  a  given 
set  of  parameters.  Although  the  above  did  not  survey  all  possible  geometric  designs  or 
flow  conditions,  it  does  demonstrate  the  ability  of  the  analysis  to  perform  this  task.  In  an 
actual  design  environment  the  experience  of  the  designer  would  be  used  to  get  a 
preliminary  design.  Then  a  series  of  cases  could  be  run  in  a  manner  similar  to  the  above, 
but  with  an  expanded  matrix  of  possibilities,  to  achieve  a  desirable  final  design.  The  use  of 
the  conjugate  heat  transfer  code  coupled  with  the  Navier-Stokes  solver  demonstrates  an 
ability  to  solve  the  more  complex  problem  where  a  surface  temperature  might  not  be 
known  beforehand.  In  many  cases  the  surface  temperature  is  not  know  or  is  not  a 
constant  value  and  can  only  be  determined  in  this  manner. 


Morton  Cases  I  -  HI 

A  series  of  three  test  cases  were  run  on  a  cylindrical  device  used  by  Morton 
International  to  deposit  SiC.  A  schematic  of  the  device  is  shown  in  Fig.  10.  Flow  enters 
from  the  top  through  the  small  injector  and  flows  into  the  large  cavity,  through  a  nominal 
throat,  impinges  onto  a  disc,  flows  over  the  disc  and  exits  through  the  bottom.  The 
injector  has  a  diameter  of  0.305  in.  and  is  15  in.  long.  The  cavity’s  diameter  is  23  in.  and 
has  a  length  of  23  %  in.  The  nominal  throat  for  the  three  cases  is  5,8  and  10  in.  and  has  a 
thickness  of  Vi  in.  The  case  shown  in  Fig.  10  is  for  the  8  in.  throat.  The  disc  always  has  a 
diameter  of  24  %  in.  and  a  thickness  of  Vi  in.  The  diameter  of  the  cavity  containing  the 
disc  is  28  %  in.  and  has  a  length  of  5  in.  The  exit  pipe  has  a  diameter  of  4  in.  with  a  length 
of  5  in.  The  overall  length  of  the  device  is  49  in.  All  dimensions  are  to  scale.  The  inlet 
injector,  the  leading  edge  of  the  first  cavity,  the  trailing  edge  of  the  second  cavity,  the  disc 
and  exit  pipe  are  assumed  to  be  adiabatic.  All  other  solid  surfaces  are  heated  to  1623  K. 
Deposition  of  SiC  also  occurs  on  all  the  heated  walls  as  well  as  the  disc.  Flow  conditions 
are  the  same  for  all  three  cases  (the  only  difference  is  the  throat  diameter)  and  are 
presented  in  Table  VII.  The  precursor  gas  is  MTS  while  the  carrier  gas  is  a  combination 
of  hydrogen  and  argon  at  the  rates  shown  in  the  table.  In  the  table  slm  stands  for  standard 
liters  per  min,  i.e.,  1000  seem.  Because  of  the  size  of  the  device  the  flow  was  assumed  to 
be  turbulent.  A  turbulent  k-s  turbulence  model  was  used.  This  resulted  in  the  solution  of 
two  additional  partial  differential  equations,  one  for  each  of  the  variables.  The 
homogenous  and  heterogeneous  chemistiy  models  are  shown  in  Table  VIII.  This 
formulation  was  obtained  from  Morton  International.  The  heterogeneous  chemistry  use  a 
global  model  where  the  rate  term  is  not  of  the  standard  Arrhenius  form. 

The  region  of  primary  interest  in  this  study  is  the  outer  radius  of  the  first,  23  in., 
cavity.  Three  computations  were  made  using  the  non-conjugate  computer  program 
corresponding  to  the  three  different  throat  diameters.  In  each  case  95  grid  points  were 
used  in  the  radial  direction  and  120  grid  points  were  used  in  the  axial  direction.  As  before 
grid  points  were  packed  in  region  where  large  gradients  of  flow  variables  were  expected. 
Streamlines  and  temperature  contours  are  shown  in  Fig.  1 1  for  the  8  in.  diameter  throat 
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case,  the  left-hand-side  showing  the  streamlines  and  the  right-hand  side  showing  the 
temperature  contours.  The  streamlines  show  two  large  recirculation  zones  in  the  first 


0.305  in 


4  in 

Figure  10.  Schematic  of  Morton  International  Device 


cavity.  The  recirculation  zone  nearest  to  the  centerline  rotates  counter-clockwise  while 
the  one  above  it  rotates  clockwise.  In  can  be  seen  that  the  recirculation  zone  nearest  the 
centerline  extends  into  the  second  cavity.  Basically  the  flow  exiting  from  the  injector 
squirts  along  the  centerline,  up  the  face  of  the  disc,  over  the  top  and  down  and  out 
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through  the  exit  pipe.  This  flow  is  only  slightly  heated  by  the  outer  walls.  In  looking  at 
the  temperature  contours  white  is  the  hottest  and  blue  the  coldest,  i.e.,  the  temperature  of 
the  fluid  entering  the  injector.  The  large  clockwise  recirculation  zone  adjacent  to  the 
outer  wall  of  the  first  cavity  is  replenished  by  this  primary  flow  with  unreacted  MTS. 
Since  flow  cannot  be  convected  across  streamlines,  this  process  occurs  by  diffusion.  The 
MTS  then  reacts  in  the  high  temperature  regions  near  the  heated  walls  and  SiC  is 
deposited  on  the  walls.  See  Fig.  1 1 . 

Fig.  12  shows  the  growth  rate  on  the  outer  wall  of  the  first  cavity  as  a  function  of 
the  throat  diameter.  The  Sin.  and  10  in.  diameter  throats  yield  growth  rates  that  are  much 
more  uniform  than  the  5  in.  diameter  throat  case.  By  comparing  the  streamlines  of  the  8 
in.  diameter  throat  in  Fig  1 1  with  the  streamlines  for  the  5  in.  diameter  throat  in  Fig  13  it  is 
easy  to  see  why  that  may  be  the  case.  While  the  8  in.  diameter  throat  case  has  two 
distinctive  recirculation  the  5  in  diameter  throat  case  has  only  one  primary  recirculation 
zone.  In  fact  the  single  recirculation  zone  is  much  more  intense  in  the  region  of  the  lower 
left  comer  of  the  first  cavity.  This  probably  leads  to  a  much  larger  growth  rate  in  this 
region  as  is  seen  in  Fig.  12.  The  second  recirculation  zone  in  the  8in  and  10  in  diameter 
throat  cases  is  moving  much  slower  than  the  single  recirculation  zone  in  the  5  in.  diameter 
throat  case  and  is  actually  going  in  the  opposite  direction.  The  direction  of  the  flow 
adjacent  to  the  surface  may  account  for  the  larger  growth  rates  for  the  8in.  and  10  in. 
cases  at  the  upper  right-hand-side  of  the  first  cavity  as  is  seen  in  Fig.  1 1  In  all  cases  the 
uniformity  of  the  growth  rate  is  less  than  that  observes  in  the  calculations  for  cases  I-IV 
above.  This  is  probably  due  to  the  size  and  design  of  the  device.  Remember  in  Cases  I-IV 
the  device  was  somewhat  sized  and  flow  conditions  modified  to  obtain  a  relative  uniform 
deposition.  No  attempt  was  made  to  do  that  in  this  case.  Instead  the  design  and  flow 
condition  was  taken  as  given  by  the  manufacturer. 


Morton  Case  I-in 


Species 

Flow  Rate  (slm) 

MTS 

7 

H2 

41 

Ar 

33 

Variable 

Value 

Poo 

200  torr 

T„ 

290  K 

Re/1« 

26572 /  cm 

He 

0.409 

Twall 

1623  K 

Table  VII.  Flow  Conditions  -  Morton  Cases 
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Homogeneous  Chemistry  -  Global  Model 

CHjSiClj  +H2  ^  ^  SiClj  +CH4  +HC1 

kb 


where  the  forward  and  backward  rates  are  given  by 

_  Ef  _ 

kf=Afe  and  kb=Abe 


where 

Af  =2.0x10^^ 
Ab  =1.1x10^^ 


cm 


mole-s 

cm 

mole^  -s 


Ef  =1.07x10^ 
Eb  =9.95x10'* 


cal 

mole 

cal 

mole 


Heterogeneous  Chemistry  -  Global  Model 


SiClj  +CH4  >H2  +2HCl+SiC(s) 


where 


and 


K  =  Aa© 


RT 


and 


where 

Af  =1.63x10^® 
=5.00x10'° 
Ab  =7.11x10® 


cm 


mole  •  s 

cm^ 

mole 

cm^ 

mole 


Ef  =7.64x10^ 
=5.16x10^ 
Eb  =7.91x10^ 


cal 

mole 

cal 

mole 

cal 

mole 


Table  VIII.  Homogeneous  and  Heterogeneous  Chemistry 
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Figure  11.  Morton  Case  8  in  Throat  Streamlines  and  Temperature  Contours 
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Axial  Location  (cm) 

Figure  12.  Morton  Case  Growth  Rate  versus  Throat  Diameter 


Figure  13.  Morton  Case  5  in  Diameter  Throat  Streamlines 


Papasouliotis  and  Sotirchos  Cases 


The  last  case  considered  was  the  experimental  device  of  Papasouliotis  and 
Sotirchos  [26].  A  schematic  of  the  device  is  shown  in  Fig.  14.  The  schematic  is  expanded 
in  the  radial  direction  by  a  factor  of  10  to  emphasize  the  substrate.  The  white  rectangle  in 
the  center  is  the  substrate.  As  with  the  other  cases  considered  the  flow  was  considered  to 
be  axisymmetric.  The  diameter  of  the  shroud/inlet  section  is  1.5  cm.  Flow  enters  from  the 
top  of  the  device  and  passes  over  the  graphite  substrate.  The  substrate  is  1.4  cm  long  and 
has  a  diameter  of  0.191  cm.  The  substrate  is  hung  from  a  digital  microbalance  that 
monitors  the  weight  of  the  substrate  as  a  function  of  time.  Overall  rates  of  deposition  are 
obtained  by  differentiating  the  weight  with  respect  to  time.  The  outer  walls  that  shroud 
the  substrate  are  heated  to  a  desired  temperature  profile.  A  region  approximately  23  cm  in 
length  is  heated  to  the  target  temperature.  Both  above  and  below  this  heating  zone  the 
temperature  rapidly  falls  off  to  the  temperature  of  the  inlet  gas.  The  temperature  was 
monitored  /measured  by  a  series  of  thermocouples  embedded  in  the  shroud  wall.  Data 
were  taken  at  a  series  of  average  shroud  temperatures,  inlet  pressures,  mass  flux  rates  and 
substrate  location.  The  substrate  location  refers  to  the  distance  of  the  midpoint  of  the 
substrate  from  the  beginning  or  top  of  the  heating  zone.  In  this  study  three  cases  were 
considered;  the  flow  conditions  are  shown  in  Table  IX.  All  cases  considered  had  a  total 
flow  rate  of  200  seem.  The  precursor  gas  is  MTS  and  the  carrier  gas  is  hydrogen.  The 
mixture  ratio  is  10%.  The  midpoint  of  the  substrate  was  always  1.8  cm  from  the  top  of 
the  heating  zone.  Calculations  were  performed  with  30  radial  grid  points  and  60  axial  grid 
points.  The  overall  length  of  the  computational  domain  was  53  cm.  The  homogeneous 
and  heterogeneous  chemistry  was  the  same  as  was  used  for  previous  MTS  cases  (see 
Tables  III  &  IV). 

I 

I  T 


53  cm 


I  T  I  t 

Figure  14.  Schematic  of  Papasouliotis  and  Sotirchos  Device 
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Comparison  of  the  computational  results  with  the  experimental  data  is  presented  in 
Table  X.  In  general  the  results  are  moderately  encouraging  at  least  with  respect  to  the 
overall  deposition  rates.  At  the  higher  temperatures  the  comparison  with  date  is  relatively 
good.  However  at  lower  temperatures  the  predicted  deposition  rate  continues  to  rise 
which  is  in  contradiction  with  the  data.  Fig.  15  shows  the  temperature  contours  for  the 
900  C  case.  Because  the  device  is  so  narrow  the  radial  direction  is  expanded  by  a  factor 
of  10.  From  this  figure  one  can  easily  see  the  heat  section  of  the  shroud.  It  is  also  easy  to 
see  that  the  temperature  (but  not  other  flow  variables)  are  essentially  one-dimensional 
above  and  below  the  substrate.  Moderate  two-dimensional  temperature  effect  can  be  seen 
in  the  region  of  the  substrate. 


Papasouliotis  and  Sotirchos  Cases  I-EQ 


Species 

Flow  Rate  (seem) 

MTS 

18.18 

H2 

181.82 

Total 

200.00 

Variable 

Value 

Poo 

100  torr 

Too 

400  K 

Re /loo 

9.41  /  cm 

Moo 

0.000398 

Twaii  (Case  I) 

900  °C 

Twaii  (Case  II) 

1000  °C 

Twaii  (Case  III) 

1075  °C 

Table  IX.  Flow  Conditions  -  Papasouliotis  and  Sotirchos  Cases 


Papasouliotis  and  Sotirchos  Cases  l-HI  Results 


Calculated  Value 

900  C 

0.04  mg/cm^-min 

0.226  mg/cm^-min 

1000  C 

0.20  mg/cm^-min 

0.131mg/cm^-min 

1075  C 

0.26  mg/cm^-min 

0. 129  mg/cm^-min 

Table  X.  Results  -  Papasouliotis  and  Sotirchos  Cases 
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Temperature  C 


0.0  300.0  600.0  900.0 

Figure  15.  Temperature  Contours  for  Papasouliotis  and  Sotirchos  900  C  case 
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•CONCLUSIONS 


A  very  general  technique  has  been  developed  to  perform  numerical  simulation  of 
CVD  processes  that  are  commonly  used  in  industry.  The  numerical  technique  has  been 
shown  to  be  very  robust  and  is  written  in  a  manner  such  that  new  physical  models  can 
easily  be  incorporated  into  the  code.  A  numerical  technique  has  been  developed  so  that 
very  low  Mach  number-high  heat  transfer  rate  cases  can  be  run  to  convergence  in 
relatively  few  numbers  of  iterations,  i.e.,  on  the  order  of  1000-2000  iterations.  Wall 
temperatures  of  an  order  of  magnitude  or  more  than  the  free  stream  temperature  can 
routinely  be  run  at  Mach  numbers  on  the  order  of  10■^  A  general  homogeneous  as  well  as 
heterogeneous  surface  chemistry  capability  has  been  developed  and  demonstrated. 

A  conjugate  heat  transfer  capability  has  been  couple  with  the  Navier- Stokes  solver 
and  results  demonstrated  for  chemically  reacting  homogeneous  as  well  as  heterogeneous 
surface  chemistry  for  low  Mach  number-high  heat  transfer  cases  with  no  appreciable 
increase  in  the  number  of  iterations  needed  to  obtain  a  converged  solution.  Compatibility 
conditions  of  temperature  and  energy  balance  are  satisfied  on  the  surface  between  the  two 
domains. 

A  series  of  six  general  classes  of  test  cases  have  been  run  and  results  generated. 
Results  have  demonstrated  that  this  technique  is  a  viable  technique  for  ‘optimizing’  a 
design  or  evaluating  an  existing  design.  When  evaluating  a  present  design  insight  can  be 
gained  about  the  physics  of  a  CVD  device  that  otherwise  could  not  be  obtained  or  if  it 
could  be  obtained  only  at  a  significant  cost.  Results  have  been  generated  and  compared 
with  experimental  data.  The  computed  results  compare  favorably  with  experimental  data 
with  respect  to  obtaining  approximate  growth  rates.  However  the  degree  of  comparison  is 
not  as  great  as  would  be  desired  and  further  work  must  be  done  in  this  area. 
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